function dz = diff_Reg2(t,z)

% Differential equations in Regime 2, where both fossil fuel and
% renewable energy are used. j>0, iF>0, iR>=0, n>=0.

global par

k = z(1);
lambda = z(2);
kB = z(3);
h = z(4);
f = z(5);
eta = z(6);
kR = z(7);
S = z(8);
N = z(9);
sigma = z(10);
nu = z(11);
qR = z(12);

dz = zeros(12,1);

%%% calculate the variables that apply to all the cases.
jkBratio = (lambda/eta/(1-par.psi))^(-1/par.psi);
jnv = jkBratio * kB;
bhdotkB = par.b * kB^(1+ par.psi) * jnv^(1-par.psi);

denom_pe = par.a*k + par.F0 - par.a * f;
pe = lambda * (par.A - par.delta) / denom_pe;

g = par.alpha0 + par.alpha1/(par.Sbar-par.alpha2/(par.alpha3+N)-S);
gdS = alpha1*(alpha3+N)^2/((Sbar-S)*(alpha3+N)-alpha2)^2;
gdN = -alpha1*alpha2/((Sbar-S)*(alpha3+N)-alpha2)^2;
gd2N = 2*alpha1*alpha2*(Sbar-S)/((Sbar-S)*(alpha3+N)-alpha2)^3;
gdSN = -2*alpha1*alpha2*(alpha3+N)/((Sbar-S)*(alpha3+N)-alpha2)^3;


c = lambda^(-1/par.gamma);
Q = par.Q0.*exp(par.popgr*t);

%%% change some denotation to simplify the expression
combo_nu = par.muR * kR * gdN;
combo_qR = par.muR*(sigma*Q-lambda*g);


if lambda == nu && lambda == qR
    
    %%% All six investments positive.
    
    LHS = zeros(5,5);
    RHS = zeros(5,1);
    
    LHS(1,:) = ones(1,5);
    LHS(2,:) = [(par.F0 - par.a*f), -(par.H0 + par.b*h), -par.a*k, ...
        -par.muR, 0];
    LHS(3,1) = (par.m+par.A-par.psi*jkBratio/(1-par.psi))*a;
    LHS(3,3) = (par.m+par.delta-par.psi*jkBratio/(1-par.psi))*a;
    LHS(4,:) = [par.a*(combo_nu+par.A-par.delta), 0, -par.a*combo_nu, ...
        denom_pe*gdN*par.muR, par.muR*kR*gd2N*denom_pe];
    
    LHS(5,:) = [combo_qR+par.A*lambda, 0, -(combo_qR+par.delta*lambda), ...
        0, -(par.A-par.delta)*k*lambda/kR];
    
    RHS(1) = par.A*k - par.m*kB - jnv - c - g*par.muR*kR;
    RHS(2) = bhotkB;
    RHS(3) = LHS(3,1)*par.delta*k + (par.A-par.delta)*par.b* ...
        jkBratio^(1-par.psi)*kB + (par.F0 + par.a*k + par.a*f) ...
        * (1/eta - 1/lambda) * jkBratio/(1-par.psi);
    RHS(4) = par.a*par.delta*combo_nu*k - denom_pe*gdSN*Q*par.muR^2* ...
        kR^2;
    
    RHS(5) = par.delta*k*(combo_qR+par.A*lambda)+ (par.A-par.delta)* ...
        sigma*Q*k*(par.muR-par.popgr*(kR*gdN)^(-1));
    
    sol = RHS\LHS;
    
    inv = sol(1);
    iB = sol(2);
    iF = sol(3);
    iR = sol(4);
    n = sol(5);
    
    
elseif lambda == nu && lambda ~= qR
    %%% iR falls to zero first. n>0.
    LHS = zeros(4,4);
    RHS = zeros(4,1);
    
    LHS(1,:) = ones(1,4);
    LHS(2,:) = [(par.F0 - par.a*f), -(par.H0 + par.b*h), -par.a*k, 0];
    LHS(3,1) = (par.m+par.A-par.psi*jkBratio/(1-par.psi))*a;
    LHS(3,3) = (par.m+par.delta-par.psi*jkBratio/(1-par.psi))*a;
    LHS(4,:) = [par.a*(combo_nu+par.A-par.delta), 0, -par.a*combo_nu, par.muR*kR*gd2N*denom_pe];
    
    RHS(1) = par.A*k - par.m*kB - jnv - c - g*par.muR*kR;
    RHS(2) = bhotkB;
    RHS(3) = LHS(3,1)*par.delta*k + (par.A-par.delta)*par.b* ...
        jkBratio^(1-par.psi)*kB + (par.F0 + par.a*k + par.a*f) ...
        * (1/eta - 1/lambda) * jkBratio/(1-par.psi);
    RHS(4) = par.a*par.delta*combo_nu*k - denom_pe*gdSN*Q*par.muR^2* ...
        kR^2;
    
    sol = RHS\LHS;
    
    inv = sol(1);
    iB = sol(2);
    iF = sol(3);
    iR = 0;
    n = sol(4);
    
elseif lambda ~= nu && lambda == qR
    %%% n falls to zero first. iR>0.
    LHS = zeros(4,4);
    RHS = zeros(4,1);
    
    LHS(1,:) = ones(1,4);
    LHS(2,:) = [(par.F0 - par.a*f), -(par.H0 + par.b*h), -par.a*k, ...
        -par.muR];
    LHS(3,1) = (par.m+par.A-par.psi*jkBratio/(1-par.psi))*a;
    LHS(3,3) = (par.m+par.delta-par.psi*jkBratio/(1-par.psi))*a;
    
    LHS(4,:) = [combo_qR+par.A*lambda, 0, -(combo_qR+par.delta*lambda),0];
    
    RHS(1) = par.A*k - par.m*kB - jnv - c - g*par.muR*kR;
    RHS(2) = bhotkB;
    RHS(3) = LHS(3,1)*par.delta*k + (par.A-par.delta)*par.b* ...
        jkBratio^(1-par.psi)*kB + (par.F0 + par.a*k + par.a*f) ...
        * (1/eta - 1/lambda) * jkBratio/(1-par.psi);
    
    RHS(4) = par.delta*k*(combo_qR+par.A*lambda)+ ...
        par.muR*(par.A-par.delta)*sigma*Q*k + ...
        par.muR*par.popgr*denom_pe*sigma*Q/par.a;
    
    sol = RHS\LHS;
    
    inv = sol(1);
    iB = sol(2);
    iF = sol(3);
    iR = sol(4);
    n = 0;
else
    %%% Fossil fuel use phases out, n = iR =0
    LHS = zeros(3,3);
    RHS = zeros(3,1);
    
    LHS(1,:) = ones(1,3);
    LHS(2,:) = [(par.F0 - par.a*f), -(par.H0 + par.b*h), -par.a*k];
    LHS(3,1) = (par.m+par.A-par.psi*jkBratio/(1-par.psi))*a;
    LHS(3,3) = (par.m+par.delta-par.psi*jkBratio/(1-par.psi))*a;
    
    RHS(1) = par.A*k - par.m*kB - jnv - c - g*par.muR*kR;
    RHS(2) = bhotkB;
    RHS(3) = LHS(3,1)*par.delta*k + (par.A-par.delta)*par.b* ...
        jkBratio^(1-par.psi)*kB + (par.F0 + par.a*k + par.a*f) ...
        * (1/eta - 1/lambda) * jkBratio/(1-par.psi);
    
    sol = RHS\LHS;
    
    inv = sol(1);
    iB = sol(2);
    iF = sol(3);
    iR = 0;
    n = 0;
    
end

dz(1) = inv - par.delta*k;
dz(2) = (par.beta + par.delta - par.A) * lambda + (par.F0-par.a*f) * pe;
dz(3) = iB - par.delta * kB;
dz(4) = kB^par.psi * jnv^(1 - par.psi);
dz(5) = iF;
dz(6) = par.beta * eta - par.b * pe * k;
dz(7) = iR-par.delta*kR;
dz(8) = Q*par.muR*kR;
dz(9) = n;
dz(10) = par.beta*sigma + lambda*par.muR*kR*gdS;
dz(11) = par.beta*nu + lambda*combo_nu;
dz(12) = (beta+delta)*qR - par.muR*(combo_qR+pe);
